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Abstract 


A new method for solving Maxwell's equations in the time domain for arbitrary values of 
permittivity, conductivity, and permeability is presented. Spatial derivatives are found by a 
Fourier transform method and time integration is performed using a second-order, semi- 
implicit procedure. Electric and magnetic fields are collocated on the same grid points, 
rather than on interleaved points, as in the finite difference time domain (FDTD) method. 
Numerical results for the propagation of a 2-D TEM mode out of a parallel plate wave guide 
and into a dielectric and conducting medium is presented. 


Introduction 


Solving Maxwell's equations by finite difference procedures in the time domain has been 
steadily gaining popularity since the pioneering work of Yee [1] and Taflove [2]; a good 
example (with more references) of this finite difference time domain (FDTD) method is 
given in the recent work of Britt [3]. Here, we present the beginnings of a new method, a 
Fourier collocation time domain (FCJLD) approach. We will concentrate here on two- 
dimensional (2-D) simulations, as these provide a good compomise between reality and 
resolution; 1-D and 3-D simulations are straightforward extensions of the results presented 
here. 

First, we will define our coordinate system and model problem, and second, set up the 2-D 
system of Maxwell's equations which we wish to solve. Third, we will explicitly discuss 
the FCTD method of solution. Fourth, we will present numerical results, and lastly, give a 
summary and conclusion. 


Problem Definition 

The physical problem we are attempting to solve by numerical simulation consists of a 
microwave horn emitting a single frequency, continuous wave signal through a dielectric 
slab and into a plasma. The motivation for this problem is the planned use of a microwave 
reflectometer to determine critical electron number densities in the stagnation region of 
hypersonic reentry vehicle. 

The 2-D model problem space is shown in Figure 1; the horizontal coordinate is z, the 
vertical coordinate is x, and the (third dimensional) y axis is perpendicular (and outward) to 
the plane of the figure. The initial conditions will consist of a sinusoidal TEM mode 
(truncated to one wavelength) in the parallel plate part of the hom. The plasma to the right 
of the dielectric slab has a gaussian profile in both the z and x directions. The physical 


dimensions of the (square) 2-D space are 10 cm on a side and the TEM mode has a 
wavelength of 2.5 cm. 


2-D Maxwell's Equations 

For now, the magnetic permeability will be assumed to have only its free space value; in 
this case, the 2-D Maxwell's equations are (SI format): 


TM eqs.: 

E d,E x = - m, 1 d z By - j x 
3jBy = 3 x Ez- 3 Z E X 
e3 t Ez= Ho^xBy'jz 


TEeqs.: 
d,B x = 3 z Ey 

£ 3 t Ey = (3 Z B x - 3 X B z ) - jy 

3 t B z = - 3 x Ey 


( 1 ) 


where the current j (and conductivity a) are: 


i: 


j(t) = a(t-x) E(x) dx ; o(t) = j- I ct(co) e i<Dt dco 




( 2 ) 


Some initial conditions which can be used for the parallel plate wave guide of Figure 1 are: 


TM modes: 

TEM: E x = By, E z =0 

TM,: ki=&k, k=^ 
2 ^ 


TE modes: 


E x = cos(k, z) cos(Pc x) 
Ez = -^ sin(ki z) sin(ik x) 


TE,: k,=&k, k = f 
B x = ^ sin(ki z) sin(^k x) 

Ey=f B x 

B z = -k cos(k, z) cos(^k x) 


( 3 ) 


In the frequency domain, the complex 'permittivity' is defined in terms of the dielectric 
permittivity e and conductivity a as: ___ 


e (co) = e + 


i q(co) 
co 


( 4 ) 


where o(co) is the Fourier transform of o(t) (and vice versa): 



a(co) e itot dco 


(5) 


but £ is not the Fourier transform of e(co) (i.e., e is always positive while Re(e(co)) need 
not be). For a plane wave of angular frequency co propagating in a plasma, a useful 
expression for o(co) is [4]: 


ofo»-iSa£i 

CO 


( 6 ) 


In the time domain, assuming a single frequency, we use for e whatever the coresponding 
value is for that frequency; for oft) in (2), we will use a constant value equal to I o(co)l at 
the frequency of interest, w 


o(t) = o 0 5(t), 00 = ^ 7 ^- 

COo 


(7) 


Thus, in (1), j= o 0 E. i n our numerical example, we will choose 


C 0 p = 2 tc x 25 GHz (maximum) 


C0o = 27txl2GHz 


so that at maximum electron density, a o - 3 mho/m. 


The FCTD Method 


The FCTD Method consists of a semi-implicit time integration scheme and a Fourier 
transform technique ("spectral method") for determining spatial derivatives. Here, spatial 
grids are NxN, where N is a power of 2; in the numerical example we will use N=128. 
Since j= c oE , the first TM equation in (1) can be written as, for example: 

£ dtE x ~ ■ Po * <J(>E X (g) 

When numerically integrating this, we can treat the last term on the right-hand-side 
implicitly, i.e., 

eEr^eES-At^feBpl + aoET 1 ] => gm-u £E ^~ At ^fe B y) 

e + At a 0 ( 9 ) 

This first-order time integration is absolutely stable with respect to the conductivity; in fact, 
any non-negative value of a is allowed on the model grid. 

Following this example, a second-order time integration mehtod for the TM equations in 
(1) is: 


Predictor 

e,e;-<M3 2 b;) 

2 

£r + 4L|o 
2 e ° 

Bf ,n =B;+&.(a 5 E„-9 zE n) 

e,ES + 4Lc2(a x B5) 

E5 +1/2 = 2 


Corrector 

e-Ej-Atc^By 1 *) 

e,+ At^ 

Co 

b;* 1 = b;+ it (d.Br in - a z E? 10 ) 

^1 + At c 2 fd x By~ lg ) 


E r + At? 2 - 


The 1/2 time step values of the fields are only used as intermediate values and are not saved 
from one time step to the next. 

The spatial derivatives in (10) are evaluated by a Fourier transform method; first, we 
transform to k-space, find the Fourier coefficients of the derivatives of the fields, and then 
transform back to x-space: 

E ( k ) = sX E (x)e- ik x => VE(x) = ££ ikE(k)e ik x 

X k (11) 

In (1 1) (and (12) below) we form a dyadic which can be specialized to a dot or cross 
product; also, x=2n(j,k)/N (0<j,k <N-1), and k=(m,n), (N/2- 1 <m,rc<N/2); j,k,m, 
and n are integers. 

The Fourier transform method is equivalent to an N 1 * 1 order finite difference method, if we 
view a generalized difference method as a convolution: 

V E(x) = V D(x-x') E(x') D(x-x') s -k sin [k-(x-x')] 

x' N k (12) 

However, instead of taking 0(N 2 ) operations for x-derivatives per point on an NxN grid , 
as the corresponding (N th order) finite difference operation would, the number of 
operations (per point) is proportional to NlogN. The ratio of accuarcy to time spent 
differentiating is thus much higher than the corresponding finite difference method. (The 
generic accuracy of spectral methods is well known [5], and will not be explicitly 
considered here.) 

An FDTD method [1,2,3], has only second order spatial differencing, which makes it 
faster, but less accurate than the FCTD method described here. Another difference is that 


in the FCTD method, both E and B are defined at the same physical grid point, while in the 
FDTD method, E and B are not defined at the same physical grid point, but rather 
'leapfrog' over one another and lie on interleaved grid points. 

Thus, in the FCTD method, E and B are collocated at all grid points. Collocation, 
however, also has another meaning in FCTD. In numerical analysis, a collocation method 
is a member of the methods of weighted residuals, and is also referred to as a 
pseudospectral method [6]. A collocation method assumes a solution in terms of known 
expansion functions and requires that the solution be exact (with respect to the function 
expansion) at the grid points; here, this means that we determine our spatial derivatives 
exactly (in terms of trigonometric functions). Since we are solving rime-dependent 
equations, we have a Fourier-collocation time-domain, or FCTD, method. (If the 
equations were linear, which they are not because of the spatial variation of the material 
constants e, 4 , and a, we would have a Fonnev-Galerkin time-domain method [6].) 

There are two kinds of boundaries in the model space. First, there are interior boundaries, 
i.e., boundaries between various media within the model space. Second, there is the 
exterior boundary of the model space itself. In the numerical example to be presented here, 
the interior boundaries are linear, and at either 0, 45, or 90 degrees to the horizontal axis. 
Material properties (such as a) are allowed to change discontinuously across the interior 0 
or 90 degree boundaries, while on the 45 degree boundaries, material properties are set to 
a value halfway between those on either side of the interface (although this is not an 
essential part of the method). The exterior boundary, because of the use of a Fourier 
transform in evaluating derivatives, is periodic; this poses few problems, however, as long 
as the bulk of the model electromagnetic energy is contained naturally within the exterior 
boundary (as occurs in the numerical example to be presented). A problem with an open 
exterior boundary will require an algorithm for absorbing energy, or expansion functions 
other than trigonometric (or both), and is not addressed here. 



Numerical Results 


As a preliminary test case, the 2-D geometry shown in Figure 1 was set up. There are three 
essentially separate regions: a parallel plate microwave "horn" on the left of the figure, a 
dielectric slab of a relative permittivity of four in the middle of the figure, and a "plasma" to 
the right of the dielectric slab. The size of the numerical grid corresponding to Figure 1 is 
128x128; the physical size of the grid is 10 cm on a side. The dielectric slab begins on 
horizontal grid point 52 and and ends on grid point 84; this corresponds to about 2-1/2 cm. 
The lower comers of the horn are at (32,48) and (48,32), while the upper comers are at 
(32,80) and (48,96). The material of the hom is assumed to be solid and have a 
conductivity of 10,000 mhos/m at the transmission frequency of 12 GHz (although any 
value can be chosen). The plasma has a conductivity of 


o(z,x) = a c exp 


2 2 2 2 


-(z-Zo) /2s z - (x-Xo) /2s 


a o =3(Qm) , z 0 =8.44cm, x o =5.00cm, s z =7.81 cm, s x = 1.56 cm 


(13) 


and is confined to the right-hand- side of the dielectric slab. 

The results are shown in Figures 2 through 7. The simulation time step size was At = 10' 14 
sec, the total number of time steps was 48,000, and the cpu time/time step was 0. 164 
seconds/At (on a Cray- 2). The total energy at the beginning of the simulation has a relative 
value of unity, and decreases as time wears on. The electromagnetic energy density, which 
is displayed in the Figures, is defined as 

T BM = i(eE 2 + t irf) (14) 

Te m can also be integrated over a selected region in physical space. This is done for the 
three disjoint regions of the model space: to the left of the dielectric slab, within the 


dielectric slab, and to the right of the dielectric slab; this regional energy is given on each 
of the Figures. 

It should be noted that there is a reflected wave in Figures 3 and 4 which barely shows up 
since its energy, which is small relative to the transmitted wave's energy, falls mostly 
below the lowest contour value in those Figures. The various reflections show up much 
better in Figures 5-7, as the primary wave has been greatly attenuated in the plasma to the 
right of the dielectric slab, and the relative energy of these reflections increases. 

Conclusion 

The results, as shown in the Figures, appear to be an accurate representation of the 
propagation of an electromagnetic wave through dielectric and conducting media (a 1-D test 
case with a linear grid of 2048 points, in which a plane wave was normally incident on a 
wide dielectric slab of relative permittivity 4, reproduced analytic predictions for reflection 
and transmission amplitudes to one part in 10 4 ). Material constants may be arbitrarily 
specified on a computational grid, so that practically any physical situation can be 
simulated. The FCTD method is of intrinsically higher accuracy than the FDTD method, as 
the FCTD method has (essentially) N th order spatial differencing, while the FDTD method 
has only 2 nd order spatial differencing. While a 2-D example was presented here, it is a 
straightforward procedure to use the FCTD method for either 1-D or 3-D problems. (We 
should also note, at this point, that finite element techniques utilizing unstructured grids 
offer an alternative which may prove to be more efficacious, in general, than either FDTD 
or spectral methods.) 

Although a Fourier transform method was used to evaluate derivatives, other function 
expansions (which may be more suitable for modelling certain boundary conditions) can be 
utilized (e.g., Chebyshev polynomials); the FCTD method is thus only one example of a 
spectral method [7]. In the future, we plan to investigate the use of alternative function 


expansions, as well as to investigate techniques for approximating absorbing boundary 
conditions and techniques for modelling the propagation of signals with broad frequency 
content. Applications which can be addressed include the simulation of microwave 
reflectometers interacting with the plasma generated in the bowshock of a hypersonic 
reentry vehicle. 
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Regional energy = 1.000 = 0.000 = 0.000 

Total energy = 1.000 at 0.00 nsec 


Figure 1. 2-D model space with initial energy distribution, along with permittivity and conductivity. 
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Regional energy = 0.696 = 0.301 = 0.000 

Total energy = 0.997 at 0.08 nsec 


Figure 2. Electromagnetic energy distribution at 0.08 nsec. 



Regional energy = 0.128 = 0.858 

Total energy = 0.987 at 0.16 nsec 


= 0.001 


Figure 3. Electromagnetic energy distribution at 0.16 nsec. 



14 



Regional energy = 0.094 = 0.785 = 0 074 

Total energy = 0.953 at 0.24 nsec 


Figure 4. Electromagnetic energy distribution at 0.24 nsec. 


Regional energy = 0.082 = 0.171 

Total energy = 0.306 at 0.32 nsec 


= 0.053 


Figure 5. Electromagnetic energy distribution at 0.32 nsec. 



Regional energy = 0.031 = 0.167 

Total energy = 0.206 at 0.40 nsec 


= 0.008 


Figure 6. Electromagnetic energy 


distribution at 0.40 nsec. 



Regional energy = 0.086 = 0.093 

Total energy = 0.183 at 0.48 nsec 


= 0.004 


Figure 7. Electromagnetic energy distribution at 0.48 nsec. 



REPORT DOCUMENTATION PAGE 


Form Approved 
OMB No. 0704-0188 


Public reporting burden for this collection of information estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, 
Gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this 
collection of information, including suggestions for reducing this burden to Washington Headquarters Services, Directorate for information Operations and Reports, 1215 Jefferson 
Davis Highway Suite 1204 Arlington, VA 22202-4302, and to the Office of Management and Budget. Paperwork Reduction Project (0704-0188), Washington, DC 20503. 


1. AGENCY USE ONLY (Leave blank) 2. REPORT DATE 


October 1991 


4. TITLE AND SUBTITLE 


A Fourier Collocation Time Domain Method for 
Numerically Solving Maxwell's Equations 


3. REPORT TYPE AND DATES COVERED 

Technical Memorandum: 


5. FUNDING NUMBERS 

WU 505-64-70-01 


6. AUTHOR(S) 

John V. Shebalin 


7. PERFORMING. ORGANIZATION NAME(S) AND ADDRESS(ES) 

NASA Langley Research Center 
Hampton, VA 23665-5225 


8. PERFORMING ORGANIZATION 
REPORT NUMBER 


9. SPONSORING /MONITORING AGENCY NAME(S) AND ADDRESS(ES) 

National Aeronautics and Space Administration 
Washington, DC 20546-0001 


10. SPONSORING /MONITORING 
AGENCY REPORT NUMBER 

NASA TM-104162 



12a. DISTRIBUTION /AVAILABILITY STATEMENT 

Unclassified - Unlimited 
Subject Category 64 


12b. DISTRIBUTION CODE 


13. ABSTRACT (Maximum 200 words) 

A new method for solving Maxwell's equations in the time domain for arbitrary values 
of permittivity, conductivity, and permeability is presented. Spatial derivatives 
are found by a Fourier transform method and time integration is performed using a 
second-order, semi-impl ici t procedure. Electric and magnetic fields are collocated 
on the same grid points, rather than on interleaved points, as in the Finite 
Difference Time Domain (FDTD) method. Numerical results for the propagation of a 
two-dimensional Transverse Electromagnetic (TEM) mode out of a parallel plate wave 
guide and into a dielectric and conducting medium is presented. 


14. SUBJECT TERMS 

computational methods, electromagnetics, Fourier analysis 


17. SECURITY CLASSIFICATION 
OF REPORT 

Unclassified 

NSN 7540-01-280-5500 


18. SECURITY CLASSIFICATION 
OF THIS PAGE 

Unclassified 


19. SECURITY CLASSIFICATION 
OF ABSTRACT 


15. NUMBER OF PAGES 

18 


16. PRICE CODE 

AO 3 


I 20. LIMITATION OF ABSTRACT . 


Standard Form 298 (Rev 2-89) 

Prescribed by AN5I Std 239-18 
29B 102 











